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We study the critical properties of a two-dimensional Ising model with competing ferromagnetic 
exchange and dipolar interactions, which models an ultra-thin magnetic film with high out-of-plane 
anisotropy in the monolayer limit. In this work we present a detailed calculation of the {5, T) phase 
diagram, 5 being the ratio between exchange and dipolar interactions intensities. We compare the 
results of both mean field approximation and Monte Carlo numerical simulations in the region of low 
values of 5, identifying the presence of a recently detected phase with nematic order in different parts 
of the phase diagram, besides the well known striped and tetragonal liquid phases. We also found 
that, in the regions of the phase diagram where Monte Carlo simulations display nematic order, 
the mean field approximation predicts hybrid solutions composed by stripes of different widths. 
Another remarkable qualitative difference between both calculations is the absence, in this region 
of the Monte Carlo phase diagram, of the temperature dependency of the equilibrium stripe width 
predicted by the mean field approximation. 
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I. INTRODUCTION 

Many ultrathin magnetic films, like e.g. Co/Cu, 
Co/Au, Fe/Cu, undergo a reorientation transition at a 
temperature Tr; for temperatures below Tr the spins 
align preferentially in a direction perpendicular to the 
film, while above Tr they align in a magnetized state 
parallel to the plane of the filroi'^'^. This reorienta- 
tion transition is due to the competition between the 
in-plane part of the dipolar interaction and the surface 
anisotropy^. Furthermore, in the range of temperatures 
where the magnetization points out of the plane, the com- 
petition between exchange and dipolar interaction causes 
the global magnetization to be effectively zero but in- 
stead striped magnetic domain patterns emerge^i^i^. In 
the limit of a monolayer film, the following dimensionless 
Ising Hamiltonian emerges as a minimal model to de- 
scribe many of the relevant physical properties of those 
materials^: 



n^-5Y.s.s, + Y.-$^ (1) 

where 5 stands for the ratio between the exchange Jo > 
and the dipolar > interactions parameters, i.e., 
5 — Ja/ Jd- The first sum runs over all pairs of near- 
est neighbor Ising spins S"; = ±1 in a square lattice and 
the second one over all distinct pairs of spins of the lattice 
(every pair is counted once); is the distance, measured 
in crystal units, between sites i and j and the energy is 
measured in units of Jd- In spite of intense theoreti- 
cal work2i£ili£i£iiaiiiii^ii^iiiii^, there are still many im- 
portant open questions regarding the critical properties 
of this model. A detailed understanding of those critical 



properties is the cornerstone of the theoretical framework 
needed to explain complex phenomena in ultrathin mag- 
netic films, such as the recently observed inverse transi- 
tion in Fe on Cu filmsfi.. 

The main difficulties to analyze the critical properties 
of this model are related to the long range character of 
the dipolar interactions, which combined with the frus- 
tration derived from the competition between interac- 
tions, adds to any theoretical approach an extra degree 
of complexity. Then, even the simplest approach, namely 
mean field approximation (MF) , leads to an infinite num- 
ber of coupled equations that, except for some particular 
situations cannot be solved exactly. For instance, in an 
early work Czech and Villain- derived an exact expres- 
sion for the MF critical temperature between the disor- 
dered and the modulated (striped) phases; however, for 
sub-critical temperatures the determination of any prop- 
erty must rely on numerical solutions of the mean field 
equations or further ansatz has to be introducediS to ob- 
tain approximated solutions of the MF equations. An ex- 
ample is the temperature dependency of the equilibrium 
stripe width; being experimentally accessible^^, reliable 
theoretical predictions of this property could be very im- 
portant to understand the basic mechanisms behind the 
complex behavior observed in these materials. Though 
MF is a powerful theoretical tool, it is known that, even 
when the corresponding equations can be solved exactly, 
neglecting fiuctuations can introduce qualitative changes 
in the critical behavior. Therefore, it is important to 
compare MF predictions with those obtained by other 
methods, in order to establish the limits of validity of 
the approximation. 

A natural way to check the mean field predictions is 
to contrast them with Monte Carlo (MC) simulations. 
However, the long range order nature of the dipolar in- 



2 



teractions makes it very difficult to simulate large sys- 
tem sizes. The ground state stripe width h is the natural 
length scale in these problems. Hence, in order to avoid 
strong finite size effects, the simulations must be carried 
out for system sizes L ^ h; this restricts the simulations 
to situations in which h is much smaller than the ex- 
perimentally observed values (typical values of the stripe 
width in Fe on Cu films, for instance, are of the orde r^^i^^ 
of 1 /xm, which corresponds roughly to h = 4000 lattice 
constants). Since the ground state value of h increases 
exponentially within 6, the values of S available for simu- 
lations are about of one order of magnitude smaller then 
the realistic values. 

In this work we carry out a detailed analysis of the 
equilibrium phase diagram of this model in the {S,T) 
space for low values of S, i.e., for < S < 4, which 
corresponds to stripe widths h < 7. By extending the 
Czech and Villain^ MF approach to the low tempera- 
ture region of the phase diagram, we performed in sec- 
tion [n] a detailed analysis of the different transition lines 
between striped states, by solving numerically the MF 
equations. Those results are compared in section lTTTl with 
MC simulations that refine previous result9^ i^°i^^i^^i^"^i^^ 
and extend them to other regions of the phase diagram. 
Our results show that, at least in the analyzed region 
of the phase diagram, several discrepancies are observed 
between both phase diagrams, which are discussed in sec- 
tion [IVl 



I J SiSj 



(2) 



where 



1 if «, j are nearest neighbors 

if i = j 
r otherwise 
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A straightforward way to derive a mean field theory 
for this Hamiltonian is the usage of the variational MF 
free energy per particlei^: 



(4) 



where N — LxL is the system size, we have taken ks = I 
and the averages are taken using the independent parti- 
cle density matrix p = Hi Pi! the one particle density 
matrices are subjected to the constraints: 



II. MEAN FIELD PHASE DIAGRAM 

The Hamiltonian ^ can be rewritten as 



E 



Pi ^ 1 SiPi = rrii. 

5'. = ±l Si = ±l 

Using the local order parameters rrii as variational pa- 
rameters we obtain the free energy functional 



[(1 + mi) In (1 + rtii) + (1 - In (1 - 



(5) 



Minimizing Eq.© respect to the order parameters rrii 
leads to the set of MF equations 
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rrii 



tanh I (3 J^^ 



1,...,7V (6) 



j{k) = Jo,e-^^^-(^^'-^"°) 

i 

= 2(5(cosfca; -I- cos/cj^) — ^ cos (fc.fy ) (8) 



Assuming periodic boundary conditions, we introduce 
the Fourier transforms 



k 



Ak.fi 



where fi is the position vector of site i, rij = r 
rn r = rht and the wave vectors k are restricted to the 

k 

first Brillouin zone. Expanding the logarithms, Eq.([5]) 
(7) can be rewritten as 
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FIG. 1: Equilibrium stripe width as a function of 5 at T = 0. 
Filled circles indicate the values of delta at which two adjacent 
striped configurations take the same energy. 



consider, for instance, a vertical striped state of width h. 
We demand the solutions of Eqs. lfTT]) to satisfy the con- 



ditions: m^^+h,y) 



-m 



Vx, y and TO(3. 



TO 



yy,y' . This restricts the harmonics in Eg.© to those 
satisfying {kx, ky) = (±(2^ + l)7r//i, 0), with I an integer 
such that 21 + 1 <h. For instance, for /i = 1 we only have 
kx — tt; for h = 2 we have k^ = ±7r/2; for /i = 3 we have 
kx = ±7r/3, tt; etc. In other words, for a stripe solution 
of width h we have h independent complex amplitudes 
TOg. In order to obtain pure real solutions we must im- 
pose mt = m_g. Replacing those conditions into Eqs.jT]) 
and (fTTj) leads to a set of h non-linear algebraic equa- 
tions for the amplitudes that can be solved numerically. 
To solve those equations we must evaluate J{kx,Q) from 
Eq. im . A suitable approximation for that function is (see 
Appendix [Aj 



which has the form of a Landau expansion. From this 
expression it is immediate that a second order transition 
between the disordered phase = Vfc, and an ordered 
phase, with non-zero order parameters, happens at the 
critical temperature^: 



Jikx,0) ~ 2S{coskx + l)-kl+27r\kx\- 
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where ({x) is the Riemann zeta-function. For the anti- 
ferromagnetic solution rni^x,y) — wo(— 1)^'*"^ we have to 
compute 



Tc = max J{k) (10) 

fe 

We calculated Tc{S) by solving Eq. l|10p numerically. 
Eqs.|[6l) can now be written as 



= tanh ■ 
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(11) 

We analyzed numerically the solutions of Eqs. lfTT]) for 
temperatures T < Tc and < (5 < 4. In particular, we 
analyzed the solutions that share the symmetries of the 
different ground states, namely, antiferromagnetic and 
striped solutions. For S < 0.425 the ground state is 
antiferromagnetio^, while for 5 > 0.425 the ground state 
is composed by stripes of width h{S). For large values of 
S we have^ h{S) ^ e^^^; for small values of S the equi- 
librium values of h can be easily evaluated numerically 
by comparing the energies of different striped configura- 
tions of increasing finite system sizes (they converge very 
quickly). h{S) is shown in FigdJ where we see that it 
attains the asymptotic exponential behavior for rather 
small values of S. 

At low but finite temperatures, the local magnetiza- 
tion inside the stripes decrease, i.e., \mi\ < 1. Let us 



J(7r,7r) = -4(5- 
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(13) 

where the last term is calculated numerically. We calcu- 
lated the MF stripe solutions for h = 1, . . . , 6 for a wide 
range of values of {S,T < Tc). To discriminate whether 
they are actually minima we analyzed the second deriva- 
tives of the free energy. For every value of h we first ana- 
lyzed the stability of the solutions, that is, for every value 
of (5 we calculated the temperature Ts{S) above which non 
trivial solutions of the above described type cease to ex- 
ist. This can be done by linearizing the corresponding set 
of equations around rh^ = and demanding the condi- 
tion of non-trivial solution, i.e., zero determinant of the 
linearized equations; this leads to the expression 



r,(5,/i) = j(J,o) 



(14) 



The stability lines are shown in FiglH together with 
Tc{S). It can be seen that for large values of S the 
stability lines accumulate near the order-disorder tran- 
sition line Tc{S), implying an increasingly large number 
of metastable states as S increases. Another remarkable 
fact, is the presence of regions near Tc{S) where no striped 



4 



solutions exists (see Figl2]). In those regions another type 
solutions appear, which are composed by parallel ferro- 
magnetic stripes of different widths. We called them 
hybrid states. The hybrid states, which we denote by 
< /i"^ /i2^ . . . /i"' >, following the notation of Selke and 
Fisher-Sfl. for the axial next-nearest-neighbor Ising model 
(ANNNI), consists in the periodic repetition of a funda- 
mental pattern composed by ni stripes of width hi (with 
opposite orientation) , followed by n2 stripes of width h2 
and so on. The regions where the hybrid states appear 
are shown in the MF phase diagram presented in FiglSl 
The boundaries between ordered phases (corresponding 
to first order transitions) where determined by compar- 
ing the free energies of the different solutions (striped and 
hybrid) using EqEJ they are shown by dashed fines in 
FigEl We found that the hybrid states appear through a 
sort of branching process near the boundary between two 
stable striped solutions as the temperature approaches 
Tc from below. For instance, the transition fine between 
the striped phases h = 1 and h = 2 ends in a triple 
point where a stable phase < 1 2 > appears between 
them; as we increase S the transition line between the 
striped phase h = 2 and the hybrid one < 1 2 > bifur- 
cates in a new triple point giving rise to the appearance 
of a < 1 2^ > phase between the < 1 2 > and the h = 2 
and so on (see inset of Figl3]). As the temperature in- 
creases more complicated hybrid states proliferate (we 
just show a few of them in Fig [3] as an example), in a 
completely analogous way as in a related model, namely, 
the three dimensional Ising model with competing short 
range ferromagnetic interactions and long range Coulomb 
interactionsSi. The MF phase diagram of that model is 
very simfiar to that of the present one, the striped states 
being replaced by lamellar ones. 

Finally, we found evidences that the proliferation of 
hybrid states also happens near the boundary between 
striped phases with larger widths (for instance, 3 and 
4), but they appear very close to Tc and the computa- 
tional effort needed to obtain an accurate estimation of 
the phase boundaries becomes very high. 



III. MONTE CARLO PHASE DIAGRAM 

Different parts of the phase diagram of this model 
were analyzed by different authors using MC numeri- 
cal simulations, for small values of 6 and small system 
sizesii^ii^'i^'iiii^. In this section we extend those results 
to other parts of the phase diagram and to larger sys- 
tem sizes (in some cases), in order to obtain a complete 
description of the small-(5 phase diagram that can be com- 
pared with the MF phase diagram. 

The MC results were obtained using heat bath dynam- 
ics on iV = L x L square lattices with periodic boundary 
conditions (Ewald sums were used to handle it2,) . We an- 
alyzed the equilibrium behavior of different quantities for 
system sizes running from L = 24 to L = 84; the max- 
imum size used for each quantity were chosen according 




FIG. 2: Stability lines Ts{S, h) for the MF striped solutions at 
finite temperatures (dashed). Full line: critical temperature 




FIG. 3: Mean field phase diagram. The numbers in the low 
temperature region indicate the equilibrium stripe width h 
of the phases. Dashed lines correspond to first order phase 
transition between ordered phases,; the full line corresponds 
to the critical temperature Tc{S). Inset: zoom of the hybrid 
region near the boundary between the striped phases h = 1 
and h = 2. 



to the associated computational effort. 

The first quantity we calculated was the order-disorder 
transition temperature as a function of 6, which we called 
Ti^^((5) (analog to TdS) in the MF case). This quantity 
was determined by means of the specific heat 

CiT) = ^ {Hf) (15) 

where (• • • ) stands for a thermal average. For some values 
of 5 we also calculated the fourth order cumulant 
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V(T) = 1 



(16) 



to characterize the order of the phase transition. 

At intermediate high temperatures (close to the order 
disorder transition and above it) this system presents a 
partially disordered phase with broken orientational or- 
der called tetragonal liquiS^'^ii^^^^^^. It is characterized 
by domains of stripes with mutually perpendicular ori- 
entations forming a kind of labyrinthine structure. At 
higher temperatures these domains collapse and the sys- 
tem crosses over continuously to a completely disordered 
phase (paramagnetic), without a sharp phase transition 
between them2,. While the existence of the tetragonal 
liquid phase has been clearly established by MC simula- 
tions, it is completely absent in the MF approximation 
(we discuss this point in section \TV\ . MC simulations 
also showed recentlj^ii that for (5 = 2 an intermediate 
phase with nematic order is present between the tetrag- 
onal liquid and the striped phase. The nematic phase is 
characterized by positional disorder and long range orien- 
tational order and is consistent with one of the two pos- 
sible scenarios predicted by a continuum approximation 
for ultrathin magnetic films in Refsi^i^^. The presence 
of the nematic phase is reflected (among several man- 
ifestations) in the appearance of two distinct maxima 
at different temperatures in the specific heat, associated 
with the stripe-nematic and the nematic-tetragonal liq- 
uid phase transitions respectively. On the other hand, 
it was shown that for 5 = 1 the specific heat present 
a unique maximum, consistent with a direct transition 
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FIG. 4: (Color on-line) Order-disorder transition temperature 
ri^' (5) and striped stability lines Ts (S) for different values of 
h (see text for details). 



the closer value of L commensurated with the modula- 
tion 2h of the corresponding ground state width). The 
results are shown by triangles joined by a continuous line 
in FiglH The order of the associated phase transition 
will be discussed in section [TVl 

We next calculated the stability of the striped phases 
in different parts of the phase diagram. It was shown 
in RefiiS that the striped phases can remain in a meta- 
stable state for values of 8 such that the equilibrium 
ground state width corresponds to a different stripe 
width. Following the same procedure as in Refj^ we 



from the tetragonal liquid to the striped phase-, sug- ^^alyzed the striped staggered magnetization 



gesting that the nematic phase is only present for some 
range of values of 8. We will call Tc the temperature 
of the high temperature peak of the specific heat, when- 
ever it presents two peaks, or the temperature of the 
unique peak if only one is present (for the system sizes 
considered and between the precision of the calculation) . 

We will call Ti^-* the temperature of the low temperature 
peak of the specific heat, when it presents two peaks. 

(2) 

While the calculation of Tc is relatively easy, the cal- 
culation of T^f^ is much more complicated and subtle, 

f 2) 

as it will be discussed below. Tc was calculated for 
different values of < (5 < 4.2 using the following simu- 
lation protocol: for each value of 8 we let first the system 
thermalize at a high enough temperature (such that it is 
in the disordered phase) during 10^ Monte Carlo Steps 
(MS); after that, we calculated the specific heat for de- 
creasing temperatures, down to a temperature well below 
the transition region. For every temperature we took the 
final configuration of the previous one and discarded the 
first 2 X 10'^ MCS for thermalization and averaged over 
10^ MCS. Every curve was averaged then over 40 inde- 
pendent runs. This calculation was performed for system 
sizes L ~ 50 for all the values of (5, in order to make the 
finite size bias comparable (for every value of 5 we choose 
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where- ^h{x) = (i — mod{x, h))/h and x, y = 1, . . . , A''. 
This quantity takes the value one in a completely or- 
dered vertical striped state of width h and zero in a 
disordered state (paramagnetic, tetragonal or nematic). 
Starting from an initially ordered vertical striped state 
at zero temperature we increase the temperature up to 
high temperatures, averaging mh at every step and using 
a similar procedure as for the calculation on Tc^\ but av- 
eraging over 10^ MCS for every temperature to diminish 
metastability effects (see discussion below) . We repeated 
this calculation for different values of S for each value of 
h; for every value of S the curve mhiT) was averaged over 
25 independent runs. Typical curves mh{T) are shown in 
FigElfor different values of h and 5. From these curves we 
estimated the stability lines Ts{5) (i.e., the temperature 
above which rrih — 0) for h — 1, ...,5 and < (5 < 4.1. 
The stability lines are shown in FigH] 

As in the MF case, we observe the existence of regions 

(2) 

below Tc where no stable striped solutions exist, at least 



6 



0.5- 




0.4 



0.3 



0.2- 



0.1 







o 5=2.1,1=0.78 
A 5=2.7,1=1.12 




X 



FIG. 5: Staggered magnetization for different S and h values: 
ft = 2, 5 = 3.4 and L = 60 (x); h ^ 3, S ^ 4.1, and L = 48 
(o); h = 5, 5 = 2.9 and L = 40 (o). 



up to 5 = 3; for values of (5 > 3 the data becomes very 
noisy (probably due to finite size effects) and the large 
error bars do not allow to identify clearly those regions for 
the present system sizes. We found that the equilibrium 
phase in those regions is a nematic one, instead of an 
hybrid state, as in the MF case (we checked the possible 
presence of hybrid states in those regions and they always 
decay after a short time into a nematic one) . The nematic 
phase is characterized by an algebraic decay of the spatial 
spin-spin correlations in one of the coordinate directions 
and an exponential decay in the other. This can be study 
through the structure factor (Fourier transform of the 
correlation function) : 



FIG. 6: Structure factor for a L = 72 system with 
S = 2.1, T = 0.78, (o) and S = 2.7, T = 1.12 (o). Continuous 
lines correspond to fittings using ||T9|| : for S = 2.1, T = 0.78, 
we obtained A = 0.033, fco = 1.39, for S = 2.7, T = 1.12, 
A = 0.0149, fco = 0.947. S{0, ky) 
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S{k) 



(18) 



where 



W e"*'^'"' . Cannas et afi^ calculated an 

approximate expression for the nematic phase structure 
factor: 



A 



A 



2Vn \ik^-ko)^ + X^ (fc^ + /co)2 + A^ 

(19) 

We ran simulations for L — 72 and different S values in 
the uncovered regions. The system was slowly heated 
from zero temperature using the same process described 
above, up to a temperature in the uncovered region, 
where we calculated S{k) by averaging over 2 x 10^ MCS. 
The typical observed behavior ofS{k) is shown in FigEl 
together with Lorentzian fittings using Eq. lfT9|) ; two typ- 
ical spin configurations in corresponding regions can be 
seen in Figl7| (compare with the results of Ref.-l^). 

Following the same steps as in the MF case, we calcu- 
lated next the transition lines between different phases 



FIG. 7: Typical nematic spin configurations for L = 72. (A) 
5 = 2.10, T = 0.78; (B) S = 2.70, T = 1.12. 



at low temperature. The transition lines were obtained 
by comparing the free energies of the striped phases of 
widths h and h + 1. The system size was chosen in these 
calculations to be a multiple of both 2h and 2{h + 1). 
In what follows we will assume that the free energy of a 
meta-stable state can be obtained by following a ther- 
modynamical path (that is, a close sequence of equili- 
brated states) from a thermodynamically stable reference 
state. To calculate the free energy of a striped phase of 
width h we first computed the internal energy per spin 
u{T, 6i) = {H) /N along a quasi-static path from an ini- 
tially low temperature To up to a working temperature 
T < Ts{6i) keeping Si constant and taking the initial spin 
configuration of a given temperature as the final config- 
uration of the previous one; the value of 6i was chosen 
well separated from the border value at zero temperature 
between the striped phases of widths h and h+1. The 
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free energy was then obtained by numerically integrating 
the thermodynamic relation 



u{(3',5,)dp', (20) 



where (3 = 1/T and I3q — 1/Tq. Once we arrive to the 
final configuration at {T,6i) we perform a second quasi- 
static path at constant temperature, by slowly changing 
S, up to a final value Sf corresponding to a striped ground 
state of width h + 1. Along this path we measure the 
average exchange energy 
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From the expression 



(21) 
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where Z is the partition function, is easy to see that 
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— = Uex{f3,5). 



(23) 



Hence, the free energy along the last path can be obtained 
by numerically integrating the equation 



UeAP,S')dS' + MP,5,) (24) 



Repeating the same procedure for the striped phase 
h + 1, but following the second path in the inverse 
sense (i.e., decreasing S from Sf down to Si) we calcu- 
lated the free energy fh+i{l3,S) at the same tempera- 
ture. The transition point St{T) is obtained from the 
equation fh{0,St) — fh+i{f3,St). We calculated the first 
order transition lines between striped phases up to ft. = 6. 
The results are shown in the MC phase diagram of FigH] 
(the transition lines between the AF and the h = 1 and 
between the h = 1 and h = 2 phases were already calcu- 
lated in Refji^; we put it here for completeness). Notice 
that all the calculated transition lines are almost inde- 
pendent of T, at variance with the MF prediction. 

We next analyzed the transition temperature Ti^' be- 
tween the nematic and the striped phases. In order to 
check whether the stability lines Ts{S) can be used to es- 
timate Tc^\ we analyzed the behavior of the specific heat 
and fourth order cumulant Eqs. (fT5| and lfT6|) around the 

(2) 

regions below Tc where no striped states exist. The 

(2) 

simulation protocol used to determine Tc is completely 
unable to detect the low temperature transition at Tc^\ 
This is because free energy barriers associated with both 



FIG. 8: (Color on-line) Monte Carlo phase diagram. The 
numbers indicate the equilibrium width of the low tempera- 
ture striped phases. The order of the different phase transi- 
tion lines are indicated in the inset (see text for details). The 
continuous lines are a guide to the eye. 



transitions for the system sizes here considered are larger 
around Ti^-* than around Tc^\ as was shown in RefJ^. 
Indeed, a rough estimation of the average times needed 
by the system to jump the free energy barrier between 
the striped and the nematic phases are of the order of the 
millions of MC, thus generating a strong meta-stability 
when the average times are of the order of 10^ MCS^^ . In 
Ref.— it was shown that an accurate estimation of ri^^ 
for S = 2 requires, for every temperature, average times of 
the order of 2 x 10* MCS. However, we verified that an av- 
erage time of 5 X 10'' is enough to determine Ti^^ between 
the error bars we are using in the present calculation (al- 
though such time scales are not enough to determine the 
height of the specific heat maximum with precision and 
therefore to allow a finite size scaling analysis). In order 
to save computational effort, we used the following pro- 
cedure for fixed values of S around the regions of interest: 

(2") 

first we ran the same simulation protocol as for Tc down 
to low temperatures and repeated it for the same param- 
eter values, but heating from a low temperature up to 
high temperatures and taking as initial configuration the 
ordered striped state. In both cases we calculated the 
internal energy u{T) along the path. This allowed us to 

determine the approximated location of Tc^\ by looking 
at the temperature range where the internal energy ex- 
hibit meta-stability^. Then we calculated C and V for 
a limited set of temperatures in that region, by taking 
averages for each temperature over a single MC run of 
5 X 10'' MCS. In order to get a more accurate estima- 
tion of Ti^' for the same values of S we also repeated the 
latter calculation for temperatures around the previous 
estimation of Tc taking averages over 10' MCS. These 
calculations were performed for S = 2.1 and 2.25 (near 
the h = 1 — h = 2 border with L = 48. The behavior of 
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FIG. 9: Specific heat and fourth order cumulant as a function 
of T for 5 = 2.25 and L = 48. 



C and V for (5 = 2.25 is shown in Figl9] (compare with 
the results of Ref i^) . We verified that the location of the 
low temperature peak of C coincides between the error 
bars with the value of Ts{S) for the same values of 6. The 
values of ri^^ for (5 = 2 (from RefM), S = 2.1, 2.25 (from 
the above calculation) and S = 2.7, 2.8, 2.9 (estimated 
from the stability lines) are shown by diamonds in the 
MC phase diagram of FigUl Reliable calculations of ri^^ 
for larger values of 6 would require system sizes that are 
out of the present computational capabilities. 



IV. DISCUSSION 

We have presented a detailed calculation of the finite 
temperature phase diagram of the Ising dipolar model in 
the range < 5 < 4, which allow striped ground state 
configurations of width up to /i = 7. We compared the 
predictions of MF approximation with extensive MC nu- 
merical simulations. Although the overall appearance of 
both phase diagrams looks similar, several differences are 
remarkable. 

The first difference to be noticed is the absence of ne- 
matic and tetragonal order in the MF approximation. 
This results from the fact that both phases are spatially 

disordered, which implies that (^S^'^ = rh^ = Wk. The 

characteristic features of those states, namely, the broken 
rotational symmetry of the nematic state and the discrete 
rotational symmetry of the tetragonal state can only be 



observed when looking at the behavior of the spatial 
correlations, or equivalently, of the structure factor— 
Eq. lfT9|) . Since fiuctuations are neglected in the MF ap- 
proximation it follows that S{k) — m^m_j: = and 
therefore the only possible spatially disordered solution 
within this approximation is the paramagnetic one. On 
the other hand, the MF approximation presents hybrid 
states solutions in the regions of the phase diagram where 
MC predicts only nematic order. Moreover, we verified 
that hybrid states are unstable in that regions, suggesting 
that (in the language of renormalization group) fluctua- 
tions play the role of a relevant scaHng field that turns the 
MF hybrid fixed points unstable towards nematic attrac- 
tors (in some sense, the hybrid states could be the closest 
state to a nematic one that can obtained when fiuctua- 
tions are neglected). This would be consistent with the 
fact that fluctuations, when included, can modify the 
continuous nature MF prediction for the phase transi- 
tion between the high temperature disordered phase and 
the low temperature ordered one: Hartree approxima- 
tion applied to the continuous version of Hamiltonian 
^ predicts a fluctuation induced flrst order transition 
any finite value of S^-, which continuously fades out for 
increasing values of 5^. In fact, MC simulations show a 
more complex scenario, where the nature of the order- 

(2) 

disorder phase transition at Tc depends on the value of 
S. 

Rastelli et al have shown that for S = the transi- 
tion is indeed continuous and belongs to the universality 
class of the nearest neighbors Ising modeU^. They also 
found a rather clear evidence of a second order transition 
for S = 0.85, but with an unusual value for the critical 
exponent f]^. However, Cannas et al have shown that 
for (5 = 1 the system presents a weak first order phase 
transitionii. These results are consistent with the pres- 
ence of a second order transition line for small values of 
S, that joins with continuous slope a first order transition 
line for larger values of (5 at a tricritical point somewhere 
between 6 — 0.85 and S — I and the unusual critical ex- 
ponents at (5 = 0.85 is probably due to a crossover effect 
near the tricritical point. There are also clear evidences 
that the transition is first order for (5 = 2 (Refs.— >i^) and 
S = 1.7,2.5 (Refpi^). The behavior of the fourth order 
cumulant observed in the present work for 6 = 2.1 and 
S — 2.25 is also consistent with a first order transition. 
For (5 = 3 the results of Rastelli et al^^ appear to suggest 
that the transition becomes continuous again. However, 
this is a matter of debate^^ and numerical results using 
a completely different technique^^ for S = 4.45 are also 
consistent with a first order transition. 

We also presented numerical evidences of the presence 
of an intermediate nematic phase between the disordered 
and the striped ones in different parts of the phase dia- 
gram. Although it seems that the nematic phase is only 
located near the border lines between striped states, the 
presence of this phase in other regions in narrow ranges 
of temperatures cannot be excluded and larger system 
sizes should be required to clarify this. 
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The existence of both type of scenarios for relatively 
large values of 6, namely one direct first order transi- 
tion from the striped phase to the tetragonal Hquid or 
two phase transitions with an intermediate nematic phase 
would be in qualitative agreement with theoretical pre- 
dictions based on a continuous approach22i2^. It is worth 
mention that Abanov et aP^ conjectured a second order 
nematic-tetragonal phase transition; since their whole 
analysis is based on mean field arguments, the disagree- 
ment with the MC results can be understood from the 
fluctuation-induced nature of the transition (see the dis- 
cussion in Refsji^ii^). Regarding the order of the tran- 
sition at Tc^\S), the situation is less clear. Cannas et 
aU^ have shown for (5 = 2 that, even the flnite size scal- 
ing is consistent with a flrst order transition, the energy 
changes continuously at Tc^\S) in the thermodynamic 
limit; this produces a saturation in the associated spe- 
ciflc heat peak, behavior that strongly resembles that 
observed in a Kosterlitz-Thouless transition. That could 
be indicative of the emergency of smectic order between 
the nematic and the striped phases for larger system sizes 
and would be in qualitative agreement with theoretical 
predictions based on a continuous approacb22iH. If that 

would be the case, probably our calculation of Tc^\S) 
overestimates the true transition temperature, since it 
is known that the speciflc heat peak locates above the 
KT transition temperature^^, and therefore the nematic 
phase would extend in larger regions of the phase dia- 
gram. However, at the present level it is very difficult to 
improve this estimation due to finite size effects. 

Regarding the low temperature behavior, a remarkable 
prediction of both MF and MC is the existence of an 
increasingly large number of striped meta-stable states 
as 5 increases. 

Finally, we found that, at variance with the MF pre- 
diction, up to (5 = 4 the transition fines between striped 
phases are completely vertical, imply temperature inde- 
pendence of the stripe width. This suggests the existence 
of some large threshold value of S, above which wall fiuc- 
tuations makes the system to cross over to a "mean field 
regime", where it starts to exhibit temperature depen- 
dency of the stripe width. 
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APPENDIX A 



When ky — Eq.|[8]) can be written as 



with 



J(fc„0) = 2,5(cosfc, + l)-5(fc,,0) (Al) 



S{kx,Q) = ^^cos{kxXi) (A2) 



r' 
i u 



where Xi is the x component of . This last term can be 
rewritten as fofiows: 



S{k^,Q) = 2'Y^ cos {k:,x)R{x) + 2C{i) (A3) 



with 



R{x) 



C(3) 



V=-oa 



(X^ -I- J/2) 3 



oo ^ 



1.202 



(A4) 



(A5) 



where x, y represents the Cartesian coordinates of each 
site and C(a:) is the Riemann zeta-function. R{x) can be 
approximated by^ 



R{x) 



dy 



2 I ,,2^ _7t2 



-cc- (.i;2 + y'^)i 



(A6) 



which has an error of %1 in the worst case (x — 1). 
Inserting this in ljA3P we get 



^(fc„0) « 4^ 



cos [kxx) 



2C(3) 



27r^ 



= fc^-27r|fc,| + — + 2C(3) (A7) 
which replaced in Eq. l|Al[) gives the expression l|12p . 
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